EXAMPLE #2: Closed Economy Sidrauski Model with Adjustment Costs. See paper for details.
Manoj Atolia: November 14, 2006.
NOTE:	 All entries concerning the determination of nonlinear solution that possibly may need to be changed are colored red. Besides these Initialization and Steady State sections are also example specific.

Off[General::spell]; Off[General::spell1]; Off[Unset::norep]; 
    Off[FindMinimum::lstol]; Off[FindMinimum::reged]; 
    Off[InterpolatingFunction::dmval]; Off[Precision::precsm];

(* 
          This cell sets the precision of calculations. *) 
{precision  =  
            100,  $MinPrecision = digits = 
              precision, $MaxExtraPrecision = Infinity};

Initialization

SetDirectory["c:"];

{τ = (3/4), ρ = (1/10), β = (1/2), δ = 
            (1/20), Ω = 2}//N

{(Subscript[bbyY, o]) = (25/100), 
        (Subscript[mbyY, o]) = (15/200), (Subscript[π, o]) = (1/2), (Subscript[θ, o]) = (1/2)}//N

{(Subscript[z, o]) = 1, (Subscript[k, o]) = 1}//N

{(Subscript[z, f]) = (6/5) (Subscript[z, o]), ω = (1/5)}//N

Steady State

U[c_, 
          m_] = (1/(1 - (1/τ))) ((c^(1 - (1/τ))) + 
              g (m^(1 - (1/τ)))); 
      Y[k_, z_] = 
        z ((((Subscript[a, 0]) + (Subscript[a, 1]) (k^(((β - 1)/β)))))^((β/(β - 1)))); ψ[
          x_] = (n/2) (((x - δ))^2); 
      n = (1/(Ω δ)); 
      (Subscript[Y, k])[k_, z_] = D[Y[k, z], k];

(Subscript[Y, o]) = 
        Y[(Subscript[k, o]), (Subscript[z, o])]; {(Subscript[c, o]) = (Subscript[Y, o]) - δ (Subscript[k, o]), (Subscript[m, o]) =  (Subscript[mbyY, o]) (Subscript[Y, o]), 
        b = (Subscript[bbyY, o]) (Subscript[Y, o]), L = (Subscript[π, o]) (Subscript[m, o]) - ρ b}; 
      g = ((Subscript[π, o]) + ρ) (((((Subscript[m, o])/(Subscript[c, o]))))^(1/τ)); {(Subscript[a, 0]), 
          (Subscript[a, 1])} =. ; {(Subscript[a, 0]), (Subscript[a, 1])} = {(Subscript[a, 0]), (Subscript[a, 1])} /. 
          FindRoot[{(Subscript[Y, k])[(Subscript[k, o]), 
                  (Subscript[z, o])]  ==  ρ + δ, (Subscript[θ, o]) == (((ρ + δ) (Subscript[k, o]))/(Subscript[Y, o]))}, {(Subscript[a, 0]), 1}, {(Subscript[a, 1]), 1}, 
            MaxIterations -> 1000, AccuracyGoal -> digits, 
            WorkingPrecision -> digits];

(Subscript[k, f]) =. ; 
      (Subscript[k, f]) = (Subscript[k, f]) /. 
          FindRoot[{(Subscript[Y, k])[(Subscript[k, f]), (Subscript[z, f])] == ρ + δ}, {(Subscript[k, f]), 
              (Subscript[k, o])}, MaxIterations -> 1000, AccuracyGoal -> digits, 
            WorkingPrecision -> digits]; (Subscript[Y, f]) = Y[(Subscript[k, f]), (Subscript[z, f])]; 
      (Subscript[c, f]) = (Subscript[Y, f]) - δ (Subscript[k, f]); 
{(Subscript[m, f]), (Subscript[π, f])} =. ; {(Subscript[m, f]), (Subscript[π, f])} = {(Subscript[m, f]), (Subscript[π, f])} /. 
          FindRoot[{g (((((Subscript[m, f])/(Subscript[c, f]))))^(-1/τ)) == (Subscript[π, f]) + ρ, L == (Subscript[π, f]) (Subscript[m, f]) - ρ b}, {(Subscript[m, f]), 
              (Subscript[m, o])}, {(Subscript[π, f]), (Subscript[π, o])}, AccuracyGoal -> digits, 
            WorkingPrecision -> digits]; 
{{(Subscript[m, o]), 
              (Subscript[c, o]), (Subscript[k, o]), (Subscript[Y, o]), b, L, (Subscript[z, o]), (Subscript[π, o]), (Subscript[c, o])/(Subscript[m, o]), 
              (Subscript[m, o])/(Subscript[Y, o])}, {(Subscript[m, f]), (Subscript[c, f]), (Subscript[k, f]), (Subscript[Y, f]), b, L, (Subscript[z, f]), (Subscript[π, f]), 
              (Subscript[c, f])/(Subscript[m, f]), (Subscript[m, f])/(Subscript[Y, f])}} // N // MatrixForm

Global Nonlinear Saddle Path

mdot = L + r b - pi m;cdot = τ c (r - ρ);kdot = i[c, k, z] - δ k;zdot = ω ((Subscript[z, f]) - z);sysvarsubs={m->m[t],c->c[t],k->k[t],z->z[t]};system={mdot,cdot,kdot,zdot};

i[c_, k_, z_]=i/.Solve[{ 
                      Y[k, z] - c - i - ψ[i/k] k == 0}, {i}, 
                    WorkingPrecision -> digits][[1, 1]];idot = 
              Dt[i[c, k, z]] /. {Dt[c] -> cdot, Dt[k] -> kdot, 
                  Dt[z] -> zdot};
eq1 := ψ''[
                      i[c, k, z]/k]/(1 + ψ'[i[c, k, z]/k]) i[
                      c, k, z]/k ((idot/i)[c, k, z] - (kdot/k)) == 
              r + δ - (((Subscript[Y, k])[k, z] + ψ'[i[c, k, z]/k] i[c, k, z]/k - ψ[i[c, k, z]/k])/(1 + ψ'[i[c, k, z]/k])); 
          eq2 := r + pi == g ((((m/c)))^(-1/τ));
{r [m_, c_, k_, z_], 
                pi[m_, c_, k_, z_]} =  {r, 
                  pi} /. Solve[{eq1, eq2}, {r, pi}, 
                    WorkingPrecision -> digits][[1]];endovarsubs={r -> r[m, c, k, z], pi -> pi[m, c, k, z]};

i[(Subscript[c, f]), (Subscript[k, f]), (Subscript[z, f])] - δ (Subscript[k, f])
(* If not zero, incorrect solution was picked for investment, i, in previous cell so change "[[1, 1]]" to "[[2, 1]]" or vice versa to pick correct solution of the quadratic equation for investment. *)

sys=Thread[{m'[t],c'[t],k'[t],z'[t]}==(-system /. endovarsubs /. 
                      sysvarsubs),List];ini:=Thread[{m[0],c[0],k[0],z[0]}=={min,cin,kin,zin},List];
  
box:= m[t]<((Min[(Subscript[m, o]),(Subscript[m, f])])/10)||m[t]>2Max[(Subscript[m, o]),(Subscript[m, f])]|| c[t]<((Min[(Subscript[c, o]),(Subscript[c, f])])/10)||c[t]>2Max[(Subscript[c, o]),(Subscript[c, f])]|| k[t]<((Min[(Subscript[k, o]),(Subscript[k, f])])/10)||k[t]>2Max[(Subscript[k, o]),(Subscript[k, f])]|| z[t]<((Min[(Subscript[z, o]),(Subscript[z, f])])/10)||z[t]>2Max[(Subscript[z, o]),(Subscript[z, f])];
  
printandplot:=(Print["Extreme Left Trajectory = ", 
                  j1, "   Closest Trajectory = ", 
                  jmin, "   Extreme Right Trajectory = ", 
                  j2,  "    No. of Subdivisions of the ϵ-ring = \
", ntraj];Print["Current Minimum Distance = ", 
                  gmindis, "  Time of Closest Approach = ", gtmin];
ParametricPlot[{{k[t],z[t]}/.gminsol,{(Subscript[k, o])+(((Subscript[k, f])-(Subscript[k, o]))/100)Sin[
                                  t],(Subscript[z, o])+(((Subscript[z, f])-(Subscript[z, o]))/100) Cos[t]}},{t, 0, gtmin},PlotRange -> All];);
  
sysdiff[t1_, t2_]:={Max[Abs[Table[((m[gtmin - t]/.gminsol)-(m[gtmin1 - t]/.gminsol1))[[1]],{t, t1, t2, 1/10}]]],Max[Abs[Table[((c[gtmin - t]/.gminsol)-(c[gtmin1 - t]/.gminsol1))[[1]],{t, t1, t2, 1/10}]]],Max[Abs[Table[((k[gtmin - t]/.gminsol)-(k[gtmin1 - t]/.gminsol1))[[1]],{t, t1, t2, 1/10}]]],Max[Abs[Table[((z[gtmin - t]/.gminsol)-(z[gtmin1 - t]/.gminsol1))[[1]],{t, t1, t2, 1/10}]]]};

distance[sol_, tsol_]:=If[firstiter == 1,
Block[{dis, dislast, tdis, steps, tbeg, tend, intsz},dis =  (10^7);dislast = (10^8);steps = 4 Floor[tsol] ;tbeg = 0;tend = tsol;intsz = tend - tbeg;
While[ (dislast - dis) > tol (10^(-5)),dislast = dis;traj=Table[Norm[Flatten[{k[t]-(Subscript[k, o]),z[t]-(Subscript[z, o])}/.sol],2],{t, tbeg, tend, (intsz/steps)}];
dis = Min[traj];tdis = 
                            tbeg + (Position[traj, dis][[1, 1]] - 
                                    1) (intsz/steps);tbeg = Max[tdis - (intsz/steps), 0];tend = Min[tdis + (intsz/steps), tsol];intsz = tend - tbeg;steps = 40];{dis, tdis}],Block[{dis, tdis},tbeg = 0;tend = tsol;min=FindMinimum[First[Norm[Flatten[{k[t]-(Subscript[k, o]),z[t]-(Subscript[z, o])}]]/.sol],{t, 9  tsol/10, tbeg, tsol},WorkingPrecision -> digits,AccuracyGoal -> 5 - Log[10, tol]];{dis, tdis} = {min[[1]], 
                        min[[2, 1, 2]]}]];

(* MESSAGES - MESSAGES. *)
diffprec := 
            Print["Difference in values of variables due to latest increase \
in precision = ",  sysdiff[0, T] // N];
diffeps := 
            Print["Difference in values of variables due to latest decrease \
in epsilon = ", sysdiff[0, T] // N];
satisprec= "THE PRECISION IS SATISFACTORY. The variables changed by less\! than\! 'tol' due to latest increase in 'precision' for current ϵ-ring.";
notsatisprec="The precision is NOT satisfactory. The variables changed by greater than\! 'tol' due to latest increase in 'precision' for current ϵ-ring.";
notsatiseps="The epsilon rings are NOT tight enough. The variables changed by greater than\! 'tol' due to latest decrease in 'epsilon'.";
increaseprec:=Print["Increased 'precision' by 4 digits to ",precision," digits."];decreaseps:=Print["Decreased 'epsilon' by a half to ",N[epsilon, precision],"."];
step2="Starting Step 2: Checking for the precision of the candidate solution.";step3="Starting Step 3: Checking for the accuracy of the candidate solution.";
congratulations= "CONGRATULATIONS! The variables changed by less\! than\! 'tol' due to latest decrease in 'epsilon'.";
done = "You are done! Hope to see you in another model. Continue with the next cell.";

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
initialize:= (tol=(10^(-5));                                                  ;
precision = 16;                                                                ;
epsilon= (10^(-2))Norm[{(Subscript[k, f])-(Subscript[k, o]),(Subscript[z, f])-(Subscript[z, o])}];                 ;);

(* MAIN CELL *)
Timing[initialize;Tmax = 4000;firstring  =  01;firstrun  =  01;
While[True,  
Label[precbegin];$MinPrecision = digits = precision;ntraj = 20;j1 = 0;j2 = ntraj;gmindis = (10^100);firstiter = 1;Label[begin];
For[j = j1,j ≤  j2,j++,theta = 2  π j/ntraj;{dkin=epsilon Cos[theta],dzin=epsilon Sin[theta]};{min,cin,kin,zin}={(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f])+dkin,(Subscript[z, f])+dzin};
If[((IntegerPart[j/2] - j/2) ≠ 
                                0 || firstiter == 1),If[firstiter == 0&&Max[Abs[{kin/pmin[[1]]-1,zin/pmin[[2]]-1}]]<(10^(-precision)),precision = precision + 4;Print["Suspecting loss in precision."];Print[increaseprec];Goto[precbegin]];
sol=NDSolve[{sys, ini},{m, c, k, z},{t, 0, Tmax},\(Method -> {"EventLocator", \ "Event" -> box, "EventAction" :> \ 
                                        Throw[
                                        Null, \ "StopIntegration"], "EventLocationMethod" -> "StepBegin"}\),MaxSteps -> 100000,WorkingPrecision -> digits];tsol = sol[[1, 1, 2, 1, 1, 2]];{dis, tdis} = distance[sol, tsol];
If[firstiter == 1, 
                              Print["Trajectory = ", 
                                j, "     Duration of Solution = ", 
                                tsol // N, "     Minimum Distance = ", 
                                dis, "     Time of Closest Approach = ", 
                                tdis]];
If[dis ≤ gmindis,gmindis = dis;gtmin = tdis;jmin = j;gminsol = sol ;pmin={kin,zin};If[gmindis < tol/3,Print["Target Tolerance = \
", tol/3 // N, "   Tolerance Achieved = ", gmindis]; Goto[end] ];];];];firstiter = 0;printandplot;{j1, j2} = 2 {jmin - 1, jmin + 1};{ntraj, jmin} = 2 {ntraj, jmin};Goto[begin];Label[end];
  
If[
                  firstring == 1  &&  firstrun == 1,  
                  T = IntegerPart[2 gtmin/3] >> T, T = << T];
Which[firstring  ==  1  &&  
                      firstrun  ==  1,{gminsol, gtmin} >> solution1; Print[step2]; 
                    precision = precision + 4; increaseprec; firstrun = 0,firstring  ==  1  &&  
                      firstrun  ==  0,{gminsol1, gtmin1} = << solution1; {gminsol, 
                        gtmin} >> solution1; diffprec; 
                    If[Max[sysdiff[0, T]] > tol, Print[notsatisprec]; 
                      precision = precision + 4; increaseprec, 
                      Print[satisprec]; Print[step3]; firstring = 0; 
                      firstrun = 1; epsilon = epsilon/2; decreaseps],firstring  ==  0   &&  
                      firstrun  ==  1,{gminsol, gtmin} >> solution2;{gminsol1, gtmin1} = << solution1;If[Max[sysdiff[0, T]] > tol,firstrun = 0; precision = precision + 4; 
                          increaseprec,diffeps;Print[congratulations];Print[done];Goto[end1] ],
firstring  ==  
                        0   &&  firstrun  ==  0,{gminsol1, gtmin1} = << solution2;{gminsol, gtmin} >> solution2;diffprec;
If[Max[sysdiff[0, T]] > tol,Print[notsatisprec]; 
                          precision = precision + 4; increaseprec,Print[satisprec];{gminsol1, gtmin1} = << solution1;diffeps;If[Max[sysdiff[0, T]] > tol,{gminsol, gtmin} >> solution1; 
                                Print[notsatiseps]; firstrun = 1; 
                                epsilon = epsilon/2; decreaseps,Print[congratulations];Print[done];Goto[end1] ]]];];Label[end1];
  
{gminsol1, 
              gtmin1} = (<< solution1)[[{1, 2}]]; {gminsol, gtmin} = (<< solution2)[[{1, 2}]]; T = << T;
For[j = T + 1,j ≤ IntegerPart[gtmin1],j++,If[Max[sysdiff[j - 1, j]] > tol,Print["The solution is accurate within tolerance 'tol' upto ",j - 1," years. Continue with next cell."];Break[]]];If[j == IntegerPart[gtmin1] + 1,Print["The solution is accurate upto ",j - 1," years. Continue with next cell."]]]

tplot = 75;
pm=Plot[{If[t < 0,(Subscript[m, o]),m[gtmin - t]/.gminsol],(Subscript[m, f])},{t, -tplot/10, tplot},PlotRange -> All,AxesLabel->{"",m},AxesOrigin -> {-tplot/10, Automatic},DisplayFunction -> Identity];
pc=Plot[{If[t < 0,(Subscript[c, o]),c[gtmin - t]/.gminsol],(Subscript[c, f])},{t, -tplot/10, tplot},PlotRange -> All,AxesLabel->{"",c},AxesOrigin -> {-tplot/10, Automatic},DisplayFunction -> Identity];
pk=Plot[{If[t < 0,(Subscript[k, o]),k[gtmin - t]/.gminsol],(Subscript[k, f])},{t, -tplot/10, tplot},PlotRange -> All,AxesLabel->{"",k},AxesOrigin -> {-tplot/10, Automatic},DisplayFunction -> Identity];
pz=Plot[{If[t < 0,(Subscript[z, o]),z[gtmin - t]/.gminsol],(Subscript[z, f])},{t, -tplot/10, tplot},PlotRange -> All,AxesLabel->{"",z},AxesOrigin -> {-tplot/10, Automatic},DisplayFunction -> Identity];
pY=Plot[{If[t < 0,(Subscript[Y, o]),Y[k[gtmin - t], z[gtmin - t]]/.gminsol],(Subscript[Y, f])},{t, -tplot/10, tplot},PlotRange -> All,AxesLabel->{"",Y},AxesOrigin -> {-tplot/10, Automatic},DisplayFunction -> Identity];
pI=Plot[{If[t < 0,δ (Subscript[k, o]),i[c[gtmin - t], k[gtmin - t], 
                              z[gtmin - t]]/.gminsol],δ (Subscript[k, f])},{t, -tplot/10, tplot},PlotRange -> All,AxesLabel -> {"", "I"},AxesOrigin -> {-tplot/10, Automatic},DisplayFunction -> Identity];
pr=Plot[{If[t < 0,ρ,r[m[gtmin - t], c[gtmin - t], 
                              k[gtmin - t], z[gtmin - t]]/.gminsol],ρ},{t, -tplot/10, tplot},PlotRange -> All,AxesOrigin->{-tplot/10, ρ},AxesLabel->{"",r},DisplayFunction -> Identity];
ppi=Plot[{If[t < 0,(Subscript[π, o]),pi[m[gtmin - t], c[gtmin - t], 
                              k[gtmin - t], z[gtmin - t]]/.gminsol],(Subscript[π, f])},{t, -tplot/10, tplot},PlotRange -> All,AxesOrigin->{-tplot/10, (Subscript[π, f])},AxesLabel->{"",π},DisplayFunction -> Identity];
pstate=ParametricPlot[{k[gtmin - t], z[gtmin - t]}/.gminsol,{t, 0, tplot},PlotRange -> All,AxesLabel -> {k, z},DisplayFunction -> Identity];
Show[GraphicsArray[{{pm,pc,pk},{pz,pY,pI},{pr,ppi,pstate}}]]

Labels={"m","c","k","z","Y","I","r","pi"};InitialSteadyState={(Subscript[m, o]), (Subscript[c, o]), (Subscript[k, o]), (Subscript[z, o]), 
              (Subscript[Y, o]), δ (Subscript[k, o]), ρ, (Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]), (Subscript[c, f]), (Subscript[k, f]), (Subscript[z, f]), 
              (Subscript[Y, f]), δ (Subscript[k, f]), ρ, (Subscript[π, f])}//N;
(Dat=Table[({tt,m,c,k,z,Y[k, z],i[c, k, z],r,pi}/.endovarsubs/.sysvarsubs/.{t ->  gtmin - t}/.tt -> t/.gminsol)[[1]],{t,0,75,1/5}])//MatrixForm;
(Dat = 
            Join[{Join[{"Initial Steady State"}, InitialSteadyState], 
                Join[{"Final Steady State"}, FinalSteadyState], 
                Join[{Time}, Labels]}, Dat]) // N // 
      MatrixForm